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Abstract 

The temporal instability of a spatially periodic, par- 
allel flow of an incompressible inviscid fluid for various 
jet velocity profiles is studied numerically using Floquet 
Analysis. The transition matrix at the end of a period 
is evaluated by direct numerical integration. For verifi- 
cation, a method based on approximating a continuous 
function by a series of step functions was used. Unsta- 
ble solutions were found only over a limited range of wave 
numbers and have a band type structure. The results ob- 
tained are analogous to the behavior observed in systems 
exhibiting complexity at the edge of order and chaos. 


I. Introduction 

The results of an investigation of the temporal instabil- 
ity of spatially periodic, parallel flow of an incompressible 
inviscid fluid for various velocity profiles are presented in 
this paper. This study is motivated by a desire to in- 
crease mixing in shear flows using multiple nozzles. 

For many years the stability of jets, shear flows, and 
wakes have been extensively studied to understand the 
growth of a small, wave-like disturbance of the basic par- 
allel inviscid flow U = U(y)i. For these unbounded flows, 
it is assumed that the flow is uniform as y -* ±oo. 

The stability characteristics of spatially periodic, par- 
allel flows of an incompressible fluid ( both inviscid and 
viscous) were studied by Beaumont 1 . Beaumont con- 
sidered periodic flows with U(y) = U(y + A) for some 
period A and all y. He considered broken line triangular 
and square velocity profiles and a continuous sinusoidal 
velocity profile^ Miles^ extended this w r ork in a study 
of the collective interaction of a compressible, periodic, 
parallel jet flow. 
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In this paper, the stability characteristics of spatially 
periodic, parallel, inviscid jet flows of an incompressible 
fluid will be studied. A sketch of the nozzle geometry 
considered is shown in Fig. 1. A broken line V-shaped 
velocity profile, a continuous Gaussian velocity profile, 
and two continuous, roughly rectangular, velocity pro- 
files, will be considered. For the steepest rectangular pro- 
file, stability will be studied for nozzle spacing- to- width 
ratios of 0.527, 0.819, 2.231, and 7.183. 

An experimental study of an array of subsonic jets 
imbedded in a square network by Villermaux and 
Hopfinger 3 and Villermaux, Gagne, and Hopfinger 4 
found the jets exhibit a self-sustained oscillation. In addi- 
tion, they found the oscillations synchronized by the fluc- 
tuating pressures gradient are strongly correlated over a 
short distance. Furthermore, measurements of the oscil- 
lations display the presence of large- wavelength travelling 
waves directed from the boundaries of the network to the 
center. 

The study discussed herein may be considered to rep- 
resent a one dimensional extrapolation to an array of 
rectangular jets separated by a distance s. The analy- 
sis depends on the jets being coherent spatially. In this 
linear temporal stability analysis of perturbations about 
the mean flow, it is the collective behavior of the periodic 
parallel jet flow that determines the behavior of the flow 
field. 

The existence of unstable solutions for the broken line 
V-shaped velocity profile and the profiles studied by 
Beaumont 1 suggested the possibility of unstable solutions 
for rectangular type profiles designed to have roughly the 
same scale. Beaumont showed that for the square bro- 
ken line velocity profile all solutions are unstable. In this 
paper, we will show that for the V-shaped velocity pro- 
file only non-growing neutral mode ( c* = 0) solutions 
exist. Stable solutions were found numerically for the 
continuous Gaussian velocity profile and two continuous, 
roughly rectangular, velocity profiles. 

II. Formulation of the problem 

The basic flow is assumed to be periodic with a maxi- 
mum value U 2 and a minimum value U\. The flow with 
velocity U{y) in the x direction is assumed to be charac- 
terized by a length scale L and a velocity scale 


1 



fi(y ) = 1 - 2 ^( 77 ), for i = 3 ... 5 


( 2 ) 


AU = U 2 -Ui. 


The dimensionless variables used herein are x = x/L, 
y = y/L, and t = tAU/L . 

Let ( U ( y ), 0) be the velocity of a steady, plane-parallel 
flow, where the x-axis is in the direction of the flow, and 

U(y) = U + —~}{y) 

where 


where 


V = A(— — 1), 

7T 


(3) 


A is the nozzle spacing parameter, y goes from 0.0 to 2n, 

1 


9s(v) = 


1 + ( sinh (iEh(iy)) 


18 


u = 


U 2 + U 1 


and / is th6 velocity profile function which varies between 

± 1 . 

The flow field is perturbed by introducing wave distur- 
bances in the velocity and pressure with amplitudes that 
are a function of y. Thus 


9*(V) = t 


and 




53 ( 17 ) = exp(-(j?) 2 ). 


(4) 


( 5 ) 


( 6 ) 


(u,v,p) = (u(y),v(y),p(y))exp[i(kx -ut)} 


where 


k = 

kL, 


QL 

U) = 

“ 


AU 

u 

Q 

k ~ 

kAU 


The profile number system used f\ = sin (27 ry) and fa = 
cos(27ry). The f\ profile case was studied by Beaumont 1 . 
The corresponding value of cPh(7])/dr} 2 is 


-1296 ei(7 )) 34 e 2 (rj) 


c 

AU 


and we define c as follows 


c _ U 
C ~ AU “ AU 2' 

By definition, A: is a real positive number that repre- 
sents the wave number in the x-direction, c is the dis- 
turbance velocity , c = c*. 4- ic*, c r is the relative phase 
velocity, and u/j = kc{/2 is the amplification rate of the 
disturbance. Prom the equations of motion, if nonlin- 
ear and viscous terms are neglected, one can obtain the 
Rayleigh equation for the y-component of the perturba- 
tion velocity as follows: 



“ JbVU/ w f ~ r 

[l + ei(77) 18 

sinh(l) 2 

= C, 

, 6126j (r]) 16 e 2 (r?) 2 


[l + ei(?7) 18 

sinh(l) 2 

r. 

36 ei (77) 18 


where 


[l + eifa) 18 ] sinh(l) s 




For <P f ±{ 7 }) / dq 2 we have 


v”(y) - v(y) 


f”(v) 

(f(y)-c) 



= 0, 


(i) 


where the primes denote differentiation with respect to 
y ' 

In this paper, f(y) is periodic such that f(y+ A) = f(y) 
where the period A used herein is the sum of a scaled noz- 
zle width, u>jv = Wh/L*, and the scaled nozzle spacing, 
s = s* /L*. In this paper the period is scaled so that 
A = s 4- iyjv = ( s * 4- w^)/L* = 27r. The velocity profile 
f(y) is not any exact solution of the Navier-Stokes equa- 
tion, but it can be considered as a simple model of some 
real periodic flow. 

The three continuous velocity profiles f{y) studied 
herein are given by 


d 2 f4iv)/dr} 2 = 


-144t? 10 ( 6 O 77 4 

[1 + T ? 6 ] 3 + [ 1 + 77 6 ] 2 ' 


For dP fz(r])/dri 2 we have 

d 2 h{r}) /dr? - 4exp(-j? 2 ) - Sr ? 2 exp(- 7 j 2 )). 


Note that 


dy 2 7 r drp ' 

The form used for the function velocity profile was 
suggested by a velocity profile function used to study the 
instability of two-dimensional wakes by Monkewitz 5 . 
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The manner in which altering the nozzle spacing pa- 
rameter, A, changes the distance between nozzles will be 
discussed next. Let xi and X 2 be the roots of f(y). The 
gap between the zero crossings is $ = X 2 - Xi • Chang- 
ing A shifts the roots of f(y). The distance between the 
nozzle center lines is the period A = 27 r. The scaled jet 
width is wjv = A - s. Results will be given in terms of 
the parameter A and a parameter which is the ratio of 
the period, A = s + wn, to the scaled jet width, w n or 
the ratio of the gap width, s, to jet width, wn* Note 


If $(y) is a fundamental matrix solution of Eq. ( 7 ) 
such that 

*(0) = I 

where I is the identity matrix, then from the Floquet- 
Bloch theorem 

#(y + 27r) = 

For y = 0, fi satisfies the matrix equation 



WN W N 


A value of 3 was used for A for some of the results pre- 
sented herein. For A = 3, the corresponding velocity pro- 
files and their first and second derivatives are shown in 
Fig. 2a and 2b. Note that the gradients in these profiles 
are steep. For the steepest velocity profile, / 5 , stability 
was studied using A = 3.0, 2.3, 1.5, andl.18 which corre- 
spond to values of A /wn = 1.527, 1.819, 3.231, and 8.183. 
Decreasing A increases X/wn and s/w n. 


III. Floquet-Bloch theory 


Since the basic flow velocity profile, /(y), is periodic, 
Eq. (1) is an example of a Floquet-Bloch problem. Differ- 
ential equation, nonlinear system, and dynamical system 
textbooks frequently discuss the mathematics of solving 
Floquet-Bloch type problems 6-16 . The method has been 
applied in solid state physics 17-21 to study the stabil- 
ity of a space vehicle 22 and a helicopter rotor 23-25 . and 
to spatially periodic flow 1 * 26-29 . A survey of the spa- 
tially periodic flow literature is presented by K. Got oh 
and M.Y. Yamada 30 . The paper by Beaumont 1 and the 
description of the Floquet-Bloch theorem by Hochstadt 9 
were particularly useful in guiding this research. 

In this section, Eq. (1) is rewritten in a form used 
to obtain the numerical solutions discussed herein. The 
second order differential equation can be described by a 
system of first order differential equations. Let 


v — Xi 
t/ = £2 

so that Eq. (1) can be rewritten as the system 


where 


X' — AX 


A<y > = ( /><») 



X is the vector X = (xi, X2) r and 


D(y) = 


' f ( y )" 
Xf { y )- c ) 



( 7 ) 


&(2n) — fx 1 = 0. 


( 8 ) 


We now introduce two solutions of Eq.(l) with initial 
values at y = 0.0 


*( 0 ) = 


<M0) 4*2 (o) 

4> i'(0) <h'( 0) 


1 0 
0 1 


( 9 ) 


Next we seek the eigenvalues, /i. of #(27r) from Eq. (8) 
as follows 


|*(2tt) - axI| = 


4>i( 2 -k)-h 4>2{ 2ir) 

4> i'(2tt) 4> 2 '{2w)- f i 


= fj . 2 - ( 4 > i { 2tt) + 4 > 2 '(2tt))^ 

+(^ > i(27t)^2 (2tt) — <fo(27r)0i ( 2n )) 
= u 2 - (<f>i (2w) + 4>2'(2n))n + 1 = 0 , 


( 10 ) 


where we have used the following relation 
M2*)4>2'{2*) -4>2(2h)4>i'(2it) = |*(2tt)| = |#(0)| = 1. 
The independent solutions of Eq. (7) have the form 

4 >( y ) = *(y)exp(^|^y) = X ( y ) exp(ry) 

The parameter T specifies the period of the eigenfunc- 
tion 4>. If r is real the eigenfunction grows or decays at 
infinity. Consequently, only imaginary values of T are 
acceptable. Thus, the eigenfunction oscillates in space 
and is called a continuous mode. The disturbance with 
Ti = 1/n, where n is a nonzero integer, has a period 
2n7r. One with Ti = 0 has the same period 27r as the 
main flow, while an irrational value of Tj means the dis- 
turbance is aperiodic. Note that the parameter T does 
not appear in the flow equation, but is attributable to 
the Floquet-Bloch theorem. 

Solutions of Eq. (7) are thus of the form 

Xi(y + 2tt) = fnX ^ y ) 

Xz(y + 2?r) = /i 2 X 2 (y) 

where the eigenvalues /ii and /i 2 represent the zeros of 
Eq. (8), provided they are distinct. 

In general, these solutions will not be periodic. 
Conditions for periodic solutions can be found by letting 
/i! = exp (i8i) and = exp (— i0t). 

Then from equation (10) 
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cos (8i) = S/2 


(ii) 

where 

<5 = 0i(2tt) + <^ / (27r) (12) 


Consequently, for a solution to be periodic, 5 must be 
real and \S\ must be smaller than 2 so that 9i is real. 

The constants fj, are termed the characteristic multipli- 
ers of the Floquet-Bloch system ( Eq.(l )) by Beaumont 1 
and the corresponding characteristic exponents are deter- 

mined by the relation r = rV + ir* = 

again in accordance with the definitions of Beaumont 1 . 

Equation (1 ) can be put in a form that shows k is an 
eigenvalue of a differential equation. 31 From Eq. 1, we 
have 


v" - k 2 


r /" 

— h i 

l(f-c)k 2 


v = 0 


(13) 


where again the primes denote differentiation with re- 
spect to y. 


Let w = v and z = v r /k and P(y, k) 




Then 


0 

1 


w 

-P 0 


z 

0 

n 


P 0 

-1 0 


0 1 


or 


Z‘ = kJH(y,h)Z 


where Z is the vector Z = (w, z) T 


J = 


0 1 
-1 0 


and 


H (y,k) 


P(v,k) 0 
0 1 


Consequently, the wavenumber k can be viewed as an 
eigenvalue of the problem. 


IV. Broken- line velocity profiles 

In this section, the stability of a broken-line V-shaped 
velocity profile and square-shaped velocity profile will 
be discussed. As shown by Beaumont 1 , one can use a 
procedure given by Rayleigh 32 to obtain exact solutions 
for various periodic broken line velocity profiles. The 
method and the nature of the solutions is demonstrated 
' using the Floquet-Bloch theory to obtain the simple ex- 
plicit solutions for a broken line profiles. Beaumont 1 dis- 
cusses additional solutions for a broken line triangular 
profile and a broken line square velocity profile. 


The Gaussian profile used herein has a resemblance to 
a V-shaped velocity profile. Consequently, in this section 
solutions for the V-shaped velocity profile are discussed. 
Since the approach is obscure more details on the solution 
are given in Appendix B. 

The rectangular profile used herein has a resemblance 
to the square velocity profile discussed by Beaumont 1 . 
The square velocity case is treated anew herein since for 
the condition of continuity of normal velocity Beaumont 1 
uses continuity of (f>/{U — c) rather than continuity of <f > . 
The condition that <f> be continuous is presented in Ap- 
pendix A where it is shown' that the normal velocity, v, 
is v = ik<j) . Both approaches yield identical results for 
the triangular velocity profile case discussed by Beau- 
mont and the V-profile treated herein which are contin- 
uous but have non continuous derivatives. However, the 
square velocity profile has discontinuities of the velocity 
profile and the stability solution presented herein is not 
the same as that derived by Beaumont 1 . Again, since 
the approach is obscure more details on the solution are 
given in Appendix C. 

For this analysis the velocity potential formulation is 
used. 33 Thus, we have 

(U - c) ( 4> " — k 2 </>) — U"<f> = 0. (14) 

In Appendix A, the flow equations in terms of the ve- 
locity potential and velocity are shown to represent two 
different views of the same problem. The conditions of 
continuity of pressure and normal velocity indicate that 
at a discontinuity of U or U' 

p(left) - p(right) (15) 

where 

P = {U- c)4t - U'<i> (16) 

and 

W(left) = W(right) (17) 

where 

V = <j> (18) 


Note again that Beaumont 1 expresses the condition of 
continuity of normal velocity using 


<t> 

(U-c) 


(19) 


The derivation of these equations is discussed in Ap- 
pendix A. 


A. Broken-line V-shaped velocity profile 


The flow profile function U ( y ) is given by 


4 


(26) 



1- %y , 0 < y < 7T 

-3+ \y , 7r < y < 2tt 


( 20 ) 


where C/(y 4- 27r) = U(y). 

Using the two sets of initial conditions at y = 0 
( (v,v') T = (1,0) T and {v,v ( ) T = (0,1) T ) ) and using 
the matching equations at y = 7r, to solve the Rayleigh 
equation (Eq.( 1)) at y = 2n yields 

4*i (2?r) = cosh(2fc7r) (21) 

4 

; cosh(&7r) sinh(fc7r) 

fc7r(l + c) 

<^ 2 , (27r) = cosh(2Ar7r) 

4 

- sinh(fc7r) cosh(fc7r). 

«7r(l + c) 

Consequently, using equation ( 12 ) and the angle sum 
relationship for sinh(a 4- /?) 

6 = 2cosh(2&7r) - (22) 


Then 


cos(0t) = cos(27iTj) 

[c 2 cosh(2fc7r) - l] 

= W^T) 

Note that solutions are available for only the neutral sta- 
bility case where c* = 0 and c = cv* For any c i n °l oqual 
to zero cos(#t) is complex and no unstable solutions exist. 

The solution given by Beaumont 1 is based on a dif- 
ferent condition for velocity continuity and this solution 
shows that the whole of the (fc, Ti) plane is unstable. 

In this section, analytic solutions were determined for 
two simple broken line velocity profiles. For these cases, 
only neutrally stable solutions were found. The numeri- 
cal procedures used to obtain solutions for the continuous 
velocity profiles will be discussed next. The continuous 
velocity profiles are scaled to resemble the simple broken 
line velocity profile. 


V. Numerical calculations 


From equation ( 11 ) 8 is real. It then follows that since 
the value of Si is zero the value of C{ is zero. Consequently, 
only neutrally stable solutions can exist. The value of 6 
for c = 0 is shown in Fig. 3a. The value of <5 is between 
±2 only over a narrow range of fc values, 0.60949283 < 
fc < 0.6573595. 

For c = 0, the value of (5 vs fc and the stability curve. 

. Ti = cos~ 1 (6/2)/2it vs fc, based on equation 11 calculated 
from the following equation 

cos(0i) = cos(27rFi) = 6/2 (23) 

2 

= cosh(2fc7r) — t — sinh(2fc7r) 

K 7T 

are shown in Fig. 3b and Fig. 3c. 

Note again that a solution only exists between fc — 
0.60949283 and fc = 0.6573595. Out side this narrow 
range of wave numbers the magnitude of <5 is greater 
than 2 and no unstable solutions to the Floquet-Bloch 
problem exist. 

B. Broken- line square shaped velocity profile 

The flow U(y) is given by: 

f 1 , 0 < y < j 

U(y) = < -1 , | <y< X (24) 

{ i , * < y < 2?r 

with U(y + 27r) = U (y) 

Using the procedure discussed in Appendix C, the fol- 
lowing equation for 6 is derived: 

s _ 2 [c 2 cosh(27rfc) - l] ^ 

c 2 — 1 


The specific point of interest is to determine for a given 
value of Ci and fc if a value of c r exists such that the 
solution to Eq. (1) is unstable. In order to determine 
the unstable regions with a numerical implementation of 
Floquet-Bloch theory one must establish a grid in the 
(ci, fc) parameter space and separately assess stability for 
each nodal point in the grid. For the results presented 
herein, the spacing in fc was 0.005 and the spacing in C{ 
WctS 0.1. 

To investigate temporal stability for a given value of c* 
and fc, iteration was used to vary the value of the phase 
velocity, Cr. At each iteration, the eigenvalues, /x, of the 
matrix ^(2?r) were then found and the corresponding 
characteristic multipliers T calculated. The iteration was 
successful if Si = 0 so that T r = 0. Also, at each iteration 
a check was made that det $(27r) = 1 . 

The iteration procedure that selects appropriate values 
of c r consisted of the following steps. First, values of c r 
in the range between —1 < c r < 1 at intervals of 0.02 
are tried. Then these results are sorted to select the two 
solutions with the smallest value of |<5i|. The two corre- 
sponding values of Cr are next used as guesses in a root 
finding IMSL subroutine (ZREAL) that tries to find two 
roots of |<5t(c r )| = 0. The solutions are then checked to 
see if |$(27 t)| = 1 and 6 r is less than 2. Both solutions 
are retained if they exist. Note that while more than two 
solutions may exist, the procedure finds at most two solu- 
tions. This procedure is expensive and time consuming, 
but alternative procedures of comparable rigor, quality, 
and generality have not been forthcoming. 

Two separate methods were used to evaluate the tran- 
sition matrix at the end of a period ( y = 27r ) so that 
6 and |^(27 t) | could be calculated. These will be dis- 
cussed next. The numerical solutions were obtained us- 
ing a CRAY YMP at the NASA Lewis Research Center. 
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A. Direct Numerical Integration 

In obtaining most of the results presented herein, equa- 
tion (1) was integrated from y = 0 to y = 27t using 
a standard Runge-Kutta procedure (IMSL math library 
routine IVPRK). Typical numerical results obtained us- 
ing numerical integration for velocity profiles /4, and 
fs using <73, 04, and 05 defined by equations 6, 5, and 4 
are given in Tables 1 through 3. Each profile is calcu- 
lated at a different pair of k and Cj values as noted in the 
Tables. 

In order to confirm the proper functioning of the com- 
puter program and assess the accuracy of the results, typ- 
ical cases were selected and verified using a spectral code 
and a method which approximated the periodic matrix 
A(y) by a series of step functions. In the next section, 
solutions obtained using the step function procedure will 
be compared to solutions obtained using direct numerical 
integration. Direct numerical integration proved faster 
than other integration procedures which were tried. In 
addition, it gives good results and copes with the steep 
gradients. 


B. Approximate Method of Determining ^(27r) 

The approximate method to be discussed consists of 
approximating the periodic matrix A(y) by a series of 
step functions. To use the method, one evaluates the 
state transition matrix by dividing a period into many 
equal parts and treats the periodic matrix A (y) as a con- 
stant matrix in each part. Consequently, the differential 
equation is treated as a constant coefficient differential 
equation in each part. The method is described by Hsu 34 
and used by Friedmann and Silverthorn to study rotor 
blade flutter 24 The period 2k is divided into K inter- 
vals denoted by k = 0, 1, 2, ...if. In the nth, interval 
the periodic coefficient matrix A (y) is replaced by a con- 
stant matrix C* = A (y K 4- ^-) where A* = y K - y K ~ x . 
Using the theory of differential equations with constant 
coefficients 35 the transition matrix of the system is given 

by 


K 

. *(2tt) = IJ exp( A k C k ) 

K = l 

The matrix exponential is defined as 

exp(A.C,) = I + f;^££ 

i=i 3 ' 

For small spatial intervals this equation converges 
quickly, and the value of the matrix exponential can be 
approximated by a sum over a finite number of terms. 
However, for this problem, it can be shown that the ma- 
trix exponential has the exact solution 35 


exp(A^C^) 


COSh^A*) 

x/ATsinh(v / ArA K ) cosh(v/^A K ) _ 


The following results are presented to provide a more 
accurate view of the calculation errors than can be con- 
veyed by the figures. 

Exact matrix exponential step function approximation 
results as a function of the number of step function in- 
tervals for the three velocity profiles studied herein are 
presented in Tables 4, 5, and 6 for the same value of k 
and Ci used in corresponding Tables 1, 2, and 3. In ad- 
dition, for velocity profile /$, the computation time in 
seconds is shown. For the Ja and fz velocity profiles re- 
sults obtained using the exact second derivative axe given 
in Tables 5a and 6a while results obtained using a numer- 
ically calculated second derivative are given in Tables 5b 
and 6b. For the step function approximation a range of 
intervals from K = 500 to K = 40, 000 were used. 

In Figures 4, 5, and 6 the discrepancy between the ap- 
proximate step function result and the direct integration 
result are presented for <$, T, and c r as a function of the 
number of intervals, K for the three velocity profiles us- 
ing the information presented in Tables 4 through 6. The 
error scales with 10 a /K b where b is approximately unity. 
For each parameter 6 , T, and Cr, the value of a increases 
as the velocity profile steepens. Figure 7 shows a plot of 
the calculation time as a function of the number of inter- 
vals, K. Note that the calculation time increases linearly 
with a slope equal approximately to 0.5. 

The integration results appear reasonable when com- 
pared to the step function approximation results. How- 
ever, using the step function approximation takes large 
amounts of computer time to get results similar to that 
obtained using integration. Consequently, the direct nu- 
merical integration procedure was used for the calcula- 
tions. Essentially, the method verified the direct numer- 
ical integration calculations. 


VI. Results 

Typical stability curves based on the numerical com- 
putation results are shown in Figs. 8-13. Results are 
presented in pairs of plots for a range of c* values. The 
upper plot shows the phase speed as a function of growth 
rate, = kd/2. The lower plot shows the periodicity 
factor T as a function of growth rate. 


A. Velocity profile effects 

In this section results obtained for the three velocity 
profiles obtained using A = 3 are discussed. For each op- 
erating condition, the unstable wave is assumed to grow 
at the maximum rate possible. While the trace of so- 
lutions are shown for a range of c* values, emphasis is 
placed on the solution with the largest growth rate. For 
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the / 3 profile (Fig. 8) solutions were found for c* in the 
range, 0.1 < d < 0.8. For the / 4 profile (Fig. 9) solutions 
were found for d in the range, 0.2 < d < 1. Similarly, 
for the /s profile (Fig. 13) solutions were found for d in 
the range, 0.1 < d < 0.92. 

For each profile the maximum growth rate occured at 
a different value of d . For the fz velocity, ( Ci) ma x — 0.4 
(Fig. 8a), for the / 4 velocity, ( Ci) m ax = 0.5 (Fig. 9a), 
and for the velocity, {d)max = 0.9 (Fig. 13a). The 
corresponding maximum growth rates are (u>i) m az = 
0.192,0.3038, and 0.2453. Thus while the range of 
( Ci)max is large, the range of (a H)max is much smaller. 
The corresponding calculated relative phase speeds are 
{c r )max = 0.1402,0.01143, and 0.2111. Thus the 
range of (cr) max is also small. 

Examining plots of the characteristic exponent, Tj, 
verses the growth rate, ( Figs 8b, 9b, and 13b) shows 
most solutions are aperiodic. Hovever, for each profile 
solutions were found with Tj near zero which have the 
same periodicity as the main flow. More importantly, 
the solution with the maximum growth rate has the same 
periodicity as the main flow. 

The trace of points in Fig 9 for d = 0.6 illustrates 
two interesting aspects of these results. First, note that 
the curve is double valued in the range 0 < < 0.12. 

Second, note that the trace consists of two islands of 
points one in the range 0 < u>i < 0.12 and the other in 
the range 2.4 < a;* < 0.28. 

Another interesting feature is the occurence of insta- 
bility waves with relative phase speeds near d — ±0.5. 
This feature can be seen in the trace of the d = 0.2 curve 
• in Fig 8 and the trace of the c» = 0.6 curve shown in Fig 
9. 


B. Effects of spacing changes 

Figures 10-13 show results obtained using velocity pro- 
file /s with four different spacings (A = 1.18,1.5,2.3, 
and 3 corresponding to the following ratios of nozzle 
spacing to nozzle width s/cj/v = 7.183,2.231,0.81 and 
0.527. Plots of the characteristic multipliers Ti and 
the phase velocity as function of growth rate for 
a range of d values are presented again as a func- 
tion of growth rate. For these spacings, the maximum 
growth rate was at the following values of ( Cj)max — 
0.3, 0.7, 0.9, and 0.9. The corresponding growth rates 
are (w<)ma* = 0.1867, 0.5547, 0.2453, and 0.2453 and 
the corresponding relative phase speeds are ( C r )max — 
0.5304, 0.1913, 0.06391, and 0.2111. The largest growth 
occurs for A = 1.5 for this velocity profile. 

The results given in Fig. 10 are unique. Results shown 
in Figs 11-13 were obtained for a large range of d values 
at least from 0.2 to 0.9. For A = 1.18 the range of d 
values was much smaller, 0.1 < d < 0.5 . Another differ- 
ence is that for the largest growth rates, for the spacing 
A = 1.5, 2.3, and 3, the relative phase velocity is in the 


range 0 < c r < 0.3 . However, for the A = 1.18 case for 
the largest growth rate c r is much larger ( c r = 0.5304). 
In addition, all solutions have positive c*-. 

The traces given in Figs 11 through 13 are similar. 
They also traces which are broken up into islands as dis- 
cussed earlier. The stability plot with A = 1.5 is inter- 
esting in that the ( d)max = 0.7 trace represented by + 
is broken up into more than four islands. 

Again in Figs 11 and 12 we note for a single trace at 
a particular value of d the occurence of instability waves 
with relative phase speeds near c r = ±0.5. These solu- 
tions switch between values' creating islands of instability 
with a particular phase speed. This is shown in Fig. 11 
for the trace c r = 0.6 and 0.7 and in Fig. 12 with the 
trace c r = 0.6, 0.7 and 0.8. 


VII. Discussion 

In this section, some general remarks about the specific 
results will be made. Then, the nature of the solutions 
will be examined. 

The collective behavior of spatially periodic, parallel, 
inviscid jet flows of an incompressible fluid supports un- 
stable waves with many common temporal stability char- 
acteristics. The results indicate that for a wide variety 
of velocity profiles and nozzle spacings the unstable wave 
with maximum growth will have a relative phase velocity 
slightly greater than zero and have the same periodicity 
as the flow. 

Extrapolation of these results to an experimental study 
of an array of subsonic jets imbedded in a square network 
by Villermaux and Hopfinger 3 and Villermaux, Gagne, 
and Hopfinger 4 suggest that the collective behavior ob- 
served might also be studied by Floquet-Bloch theory. 

Systems of equations with coefficients periodic in time 
arise in the study of the stability of helicopter rotor 
blades flapping motions 23-25 , fluid column stability in 
the presence of periodic accelerations 36 , the instability 
of columns under periodic loading 37 , the stability of os- 
cillatory Stokes layers 38 and the instability of any system 
undergoing parametric excitation. Systems of equations 
periodic in space arise in solid state physics 17-21 , spa- 
tially periodic hydrodynamic flow 1,26-29 . The question 
of interest in many cases is 

• determining the frequency of the instability which 
has the highest growth for systems undergoing pe- 
riodic excitation 

• the wavelength or the wave number of the instabil- 
ity which has the highest growth for systems with 
spatially periodic coefficients. 

Theoretical models for systems where the variation is 
sinusoidal exhibit a continuous family of stable states. 
However, for the spatial variation considered herein one 
finds in addition islands of instability. The eigenfunc- 
tions of interest are the ones that do not vanish at infinity 
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which are associated with the continuous eigenvalues and 
are called continuous modes. The new results obtained 
with steep velocity profiles indicate that these solutions 
can exist in discrete regions. 

Beaumont 1 , showed that for sinusoidal velocity profiles 
the result is a region of instability and for a square veloc- 
ity profile ail points are unstable. However, it is shown 
herein that for velocity profiles that are almost square, 
Islands of instability are found. This simple formulation 
for almost square velocity profiles produces results that 
are surprising, complicated and essentially unpredictable. 

These results can be categorized as follows. For the 
v , triangular, and sinusoidal velocity profiles, the (c*,fc) 
space is divided into two regions one stable and one un- 
stable and the space is highly ordered . For the square 
velocity profile all solutions are unstable and the space 
is highly disordered . In this paper, it is shown that by 
selecting certain velocity profiles, one observes a phase 
transition between highly ordered and completely disor- 
dered ( Ciyk ) phase space, analogous to the phase transi- 
tion between the solid and fluid states of matter. This 
analogy is also similar to the one made by Langton in 
the study of cellular automata 39,40 . Consequently, these 
results represent another example of a system exhibiting 
complexity at the edge of order and chaos 40,41 . 

The procedure used to find solutions has some similar- 
ity to that used to study the Julia sets and the Mandel- 
brot sets which are fractal geometric objects. However, 
the search for solutions at each (c*, k) point is more dif- 
ficult and the comprehensive method used requires more 
computer time. Consequently, enough work has not been 
done at this time to discern if the solutions are located in 
rigid domains or if the domain boundaries have a fractal 
nature. However, the basic procedure used involves find- 
ing the roots of |£t(cr)| = 0. Newton's method provides 
a means to compute them using the following dynamical 
system 



Discussions of the application of Julia set theory to 
Newton’s Method 42 and a brief discussion by Shub on 
the theory of the complexity of equation-solving 43 sug- 
gest the solutions might have a fractal nature. The sensi- 
tivity of the problem to spacing and the complexity of the 
results suggest that this is might be an interesting prob- 
lem for someone working on the theory of the complexity 
of equation-solving. 


VIII. Conclusions 

For a range of amplification rates, solutions for the 
. velocity profiles considered were found. For these cases, 
solutions exist over a limited range of wave numbers and 
have a band type structure. 

These results represent another example of a system 
exhibiting complexity at the edge of order and chaos 


Appendix A: Formulation of Flow Equations 

The analytical study of the broken line V-shaped ve- 
locity profile uses a velocity potential formulation of 
Rayleigh’s stability equation. The numerical study uses 
a velocity perturbation formulation of Rayleigh’s stabil- 
ity equation. In this appendix, it is shown how these 
two formulations are related. The derivations starts with 
the two dimensional equations presented by Drazin and 
Reid 33 (Eq. 21.11) written in terms of u and v rather 
than u and w since in this paper we consider U ( y ) in- 
stead of U (z ) . 

Thus we have 


ik(U — c)u 4- U'v = — ikp A.l 

ik(U - c)v = —Dp A. 2 


iku + Dv = 0 A.3 


where / and D represent the d/dy operator. Introducing 
a stream function ip such that the two components of the 
disturbance velocity are given by 



A. 4 


v = 


dip 

dx 


A. 5 


Let 


ip = <p{y) exp l *( x ct ) A.6 


Then 


u = <p' A.7 

and the normal velocity is given by 
v = —ikp A.8 

Consequently, equation A.3 is automatically satisfied. 

Substituting Equations A.7 and A.8 into Equation A.l 
and solving for p yields 


p = U’p - (17 - c)4>*. A. 9 


Substituting equation A.8 and equation A. 9 into equation 
A.2 yields a velocity potential formulation of Rayleigh’s 
stability equation 




U n 


(U-c) 


+ Jfc 2 


= 0 . 


A. 10 


The velocity formulation of Rayleigh’s stability equation 
is obtained by first solving equation A.2 for v 
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A. 11 


Using the flow profile function given by equation 20 we 
have 


v = 


Dp 

‘ ik(U-c)‘ 


Then equation A. 3 is solved for u 

= -1 

B.7 

iku — —Dv A. 12 

N | 

1 

II 

B.8 

Next one substitutes equation A. 3 into equation A.l 


U{ 7T+) = -1 

B.9 

_([/ _ c )Dv 4- U'v = —ikp A. 13 

,, , 2 

B.10 

U (tt + ) = - 

Then one takes the derivative of equation A. 13 with re- 

7 r 


spect to y which yields 

Evaluating equations B.3 thru B.6 at 

7r we have 

{U - c)(-D 2 v) 4- U u v = - ikDp A.14 

0i_ = cosh(fc7r) 

B.ll 

Last, one uses equation A. 2 to substitute for Dp to obtain 

0i _ = &sinh(/c7r) 

B.12 


D 2 v — v 


U" 

(1 U-c ) 



= 0 


A. 15 


The conditions for the continuity of pressure and nor- 
mal velocity at a discontinuity of U or U ' for the broken 
line profile given by equations 15 thru 18 are easily de- 
rived from equations A. 8 and A. 9. 


0 1+ = cosh(fc7r) B.13 

0' 1+ = fcsinh(fc7r) + M\ c k B.14 


0 2 _ = ^sinh(fc7r) B.15 

fC 

(j> 2 _ = cosh(fc7r) B.16 


Appendix B: Stability analysis for the broken line 
V-shaped velocity profile 

The analytical study of the broken line V-shaped pro- 
file uses a velocity potential formulation of Rayleigh’s 
stability equation. For the broken line profile U n is zero 
and the resulting differential equation 

0" - 0 [fc 2 ] = 0. B.l 

has solutions which are linear combinations of sinh(Ary) 
and cosh(Ary). 

We now introduce two solutions of Eq. (B.l) 0i (y) and 
0 2 (y) with initial values at y = 0.0 given by equation 9. 
In the region 0 < y < n we have solutions 

(f>i{y) = cosh(fcy) B.3 

and 

My) = - sinh(fcy). B.4 

In the region 7 r < y < 2ir we have solutions 

<f>i = cosh(fcy) 4- Mic sinh(A;(y - 7r)) B.5 

and 

02 — ~ sinh {ky) 4- M 2s sinh(A;(y - 7r)) B.6 

rC 


0 2 + = t sinh(fc7r) 

K 


B.17 

02+ = cosh(fc7r) 4- M 2s k B.18 


Equations for continuity of the normal velocity ( equa- 
tions 17 and 18) are satisfied at the velocity gradient 
discontinuity position 7r . 

Equations for continuity of pressure ( equations 15 and 
16) at the velocity gradient discontinuity position n yield 
values for M\ c and M 2s - 


4 cosh(fc7r) 
7 r (1 4- c) 


and 


M 2s = - 


4 sinh(/:7r) 
(l + cj“ 


B.19 


B.20 


At 27r we have 

0 x+(2tt) = cosh(27rfc) 4- Mi c sinh(7rfc) B.21 

and 

02+(2tt ) = cosh(27rfc) 4- M 2s kcos /i(7rfc) B.22 

Using equation 12 and the angle sum relationship for 
sinh(a 4- /?) we can obtain equation 22. 


. Next the conditions for continuity of pressure and normal 
velocity given by equations 15 thru 18 are applied. 


9 



Appendix C: Stability analysis for the broken line 
square velocity profile 


C.13 

C.14 


4> 1 , 2 + = N lc 

0i,2+ = kMic 


Again, as with the broken line V-shaped profile, the 
analytical study of the broken line square velocity profile 
uses a velocity potential formulation of Rayleigh's stabil- 
ity equation. Additionally, again, the broken line profile 
U n is zero and the resulting differential equation (B.l) 
has solutions which are linear combinations of sinh (ky) 
and cosh (ky). 

As before two solutions of Eq. (B.l) <j>\(y) and 0 2 (?/) 
axe introduced with initial values at y = 0.0 given by 
equation 9. 

In region one 0 < y < f we have solutions 
*i,i(j/) = cosh (ky) C.l 

and 

4>2,i{y) = Isinh(fcy). C.2 
In the region j < y < we have solutions 
*i ,2 = IV lc cosh(fc(y - |) + M lc sinh(k(y - |)) C.3 

and 

* 2,2 = sinh(k(y — ^)) + M\ s cosh(fc(y — ^)) C.4 

In the region ^ < y < 2tt we have solutions 

* 1,3 = N 2c cosh (k(y - — ) + M 2c sinh(fc(y - — )) C.5 

and 

* 2,3 = -iV 2 »sinh(A:(j/— — ))+M 2 s cosh (k(y~)) C.6 

Next the conditions for continuity of pressure and nor- 
mal velocity given by equations 15 thru 18 are applied. 
Note that U 9 = 0. 

Using the flow profile function given by equation 20 we 
have 



C.7 


C.8 

^ + )— 1 

C.9 

mm 

+ 

II 

o 

C.10 


Evaluating equations C.l thru C.4 at J we have 
0i t i_ = cosh(fc^-) C.ll 
<f>[ i_ = fcsinh(fc^-) C.12 


02,1- = ^sinh(fc^) C.15 
02,1- = cosh(fc^) C.16 


02+ — Mi s 

C.17 

02+ = Nls 

C.18 

Equations for continuity of the normal velocity ( equa- 
tions 17 and 18) are satisfied at the velocity discontinuity 

position j if 


Nu = cosh (*I) 

C.19 

= **<*») 
k 

C.20 


Equations for continuity of pressure ( equations 15 and 
16 ) at the velocity discontinuity position 7 r yield values 
for Mic and N\ s . 

M le = Uinh(lc^) C.21 

c+1 v 2' 

and 

N ls = cosh(*|) C.22 

Equations for continuity of the normal velocity ( equa- 
tions 17 and 18 ) and pressure ( equations 15 and 16 ) 
yield at the velocity discontinuity position ~ : 

[iV lc cosh(fc7r) + Mi c sinh(A;7r)] 

-N 2c = 0 C.23 

(—1 — c)[A^i c sinh(fc7r) + M\ c cosh(A;7r)] 

-(1 - c)M 2c = 0 C.24 


[^^sinh(Ar7r) + Mi s cosh(A:7r)] 
k 

-M 2s = 0 C.25 

(— 1 — c)[IVi, cosh(fc7r) + Mi, sinh(fe7r)] 
-(1 - c)M 2s = 0 C.26 


Solving for N 2c , M 2c , M 2s , and N 2s and using the angle 
sum relationships for cosh(a+/?) and sinh(a+/3) we find: 


N 2c 


M 2c 


n 2s 


m 2 . 


1 

1 + C 

, ,3 kn. , .kn ' 

c cosh(-^—) -I- cosh(— ) 

1 

c — 1 

. _ ,3 fc7T, . , .kit' 

csmh(— ) + smh(— ) 

1 

c — 1 

ccosh(^) - cosh(y) 

1 


csinh(^) -sinh(^ 

k{c+ 1) | 


C.27 

C.28 

C.29 

C.30 
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At 27 r we have 

7rk 

<j>\+( 2n) = N 2c cosh(— ) + M 2c sinh(— ) C.31 

and 

4>I> + (27r) = iV 2s cosh(^) + M 2s k sinh( C.32 

Using equation 12 and the angle sum relationships we 
can obtain equation 25. 
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Table 1 Integration results for (sinhfa / sinh(l))) is velocity profile 


fir 

r Cr 

e 1.5430059290031E-2 

0.2487720990379 0.9002607287015 

= 0.4, k = 0.115 


Table 2 Integration results for rf velocity profile 


<5 r 

r 

c r 

-1.795749441083, 

7.2555284510676E- 2 

2.684018292084E-2 


a = 0.5, k = 1.195 


Table 3 Integration results for Gaussian velocity profile 


fir 

r 

Cr 

-1.76193559911 

7. 84463124484582? - 2 

0.1914746533746 


Ci = 0.4, k = 0.9 


Table 4 Exact Matrix Exponential Step Function Approximation 
Results for (sinh{r] / sinh( l))) 18 velocity profile 


Intervals,^ 

fir 

r 

C r 

Time 

12000 

0.1324704157898 

0.2394506125696 

0.8954450299652 

6106 

13000 

0.1233419809349 

0.2401785210625 

0.8958327744015 

6617 

14000 

0.1155161824543 

0.2408023920132 

0.8961633816575 

7127 

15000 

0.108736216033 

0.2413427747509 

0.8964479072959 

7639 

16000 

0.1028055481471 

0.2418153838495 

0.8966953578756 

8149 

17000 

9.75739887189692? - 2 

0.24223222159 

0.8969125382767 

8658 

18000 

9.29248002178572? - 2 

0.2426026126168 

0.8971046810721 

9157 

19000 

8. 8765848422 10 IE — 2 

0.2429339134897 

0.8972758834688 

9677 

20000 

8.50234849469062? - 2 

0.2432320028322 

0.8974293902613 

10188 

30000 

6 . 1 3366923986492? - 2 

0.2451182120209 

0.898389137051 

15271 

35000 

5.4573663580242? - 2 

0.2456566231031 

0.8986594651123 

17813 

40000 

4.9502716977565 E - 2 

0.2460602929786 

0.8988610967745 

20377 

= 0.4, A: = 0.115 
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Table 5a Exact Matrix Exponential Step Function Approximation 
Results for rj 6 velocity profile using exact /" 


Intervals,!^ 

Sr 

r 

c r 

5000 

-1.731639624708 

8.3398757988445E- 2 

2.663504066365E - 2 

10000 

-1.763593751114 

7.8166992537224E — 2 

2.6737302451565J5 — 2 

20000 

-1.779569996798 

7.54271 10978254E— 2 

2.6788433789913E - 2 

21000 

-1.780330756172 

7.5294351899531E- 2 

2.6790868618323E-2 


Ci = 0.5, fc = 1.195 

Table 5b Exact Matrix Exponential Step Function Approximation 
Results for rf velocity profile using numerical evaluation of f" 


Intervals,^ 

Sr 

r 

C r 

1000 

-1.476622305991 

0.1178118119098 

2.5823569373648E — 2 

10000 

-1.764242903764 

7.8057391206773J5— 2 

2.6743768291609E - 2 

20000 

-1.78021716754 

7.5314187830809E-2 

2 .67948951 55659E - 2 

30000 

-1.785541787533 

7.4379094561713J5 - 2 

2.681 1937462023E - 2 


a = 0.5, k = 1.195 


Table 6a Exact Matrix Exponential Step Function Approximation 
Results for a Gaussian velocity profile using exact value of f" 


Intervals,^ 

<5r 

r 

C r 

5000 

-1.743858151466 

8.1434477636742E-2 

0.1914012580744 

10000 

-1.752887165675 

7.9954766508094E — 2 

0.1914378621255 

20000 

-1.757396473618 

7.9206314493282E- 2 

0.1914561345435 

25000 

-1.758297914811 

7.90559160312E-2 

0.1914597866106 

30000 

-1.758898796215 

7.8955517837813E-2 

0.1914622208468 


Ci = 0.4, k = 0.9 

Table 6b Exact Matrix Exponential Step Function Approximation 
Results for a Gaussian velocity profile using numerical evaluation of 

f" 


Intervals,/^ 

Sr 

r 

Cr 


-1.671129084883 

9.2570372553067E — 2 

0.1911055992358 


-1.743886073843 

8.1429939437306E — 2 

0.1914016096238 


-1.752915022464 

7.9950162541102E — 2 

0.1914382134012 


-1.757424297442 

7.9201676217921E — 2 

0.1914564856811 


-1.758926608932 

7.8950867920244E - 2 

0.1914625719339 


a = 0.4, k = 0.9 
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Fig. 1 Nozzle configuration. 
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a) Velocity profiles. 



b) Second derivative of velocity profiles. 


Fig. 2 Velocity profiles and second derivative of velocity profiles. 
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a) 6 versus k with c = 0. 




Fig. 3 Stability curves for broken line v shaped 
velocity profile. 
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Fig. 4 Error in as a function of the number of intervals. 
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Fig. 8 Stability plots for velocity profile withA = 3. 
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Fig. 6 Error in f ■ as a function of the number of intervals. 
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Fig. 7 Calculation time in seconds as a function of the number 
of intervals. 
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Fig. 9 Stability plots for J 4 velocity profile withA = 3. 
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